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We use the truncated Wigner approximation to derive stochastic classical field equations for 
the description of polariton condensates. Our equations are shown to reduce to the Boltzmann 
equation in the limit of low polariton density. Monte Carlo simulations are performed to analyze 
the momentum distribution and the first and second order coherence when the particle density is 
varied across the condensation threshold. 



I. INTRODUCTION 

Condensates of microcavity polaritons [l[ are a solid 
state realization of the two dimensional Bose gas. Their 
succesful creation relies on the peculiar nature of the mi- 
crocavity polariton quasiparticle, that combine a very 
light effective mass (high quantum degeneracy tempera- 
ture) with intcrparticle interactions that provide efficient 
relaxation. The formation of spontaneous coherence in 
these systems is now routinely achieved in several labo- 
ratories @, i l, i] . 

One crucial difference between polariton condensates 
and other realizations of the two dimensional Bose gas 
such as liquid 4 He films [6[ and tightly confined ultra- 
cold atomic gases 0] comes from the finite life time of 
the microcavity polaritons of the order of a few ps. In 
order to compensate for the polariton losses, new par- 
ticles can be continuously injected into the microcavity. 
The resulting steady state is not a thermal equilibrium 
one, still it shows features expected for an equilibrium 
BEC. For example, the tail of the momentum distribution 
can in many cases be fitted by an exponential Maxwell- 
Boltzmann decay. The lack of full thcrmalization is al- 
ready clear from the fact that the extracted temperature 
is in general not equal to the temperature of the reservoir 
constituted by the semiconductor lattice @, 0, IH . 

Effects that have no counterpart at equilibrium have 
been observed in polariton condensates. For example, the 
condensate state can depend dramatically on the size of 
the excitation spot [H |||: in the case of a large pump 
spot the usual condensation around zero momentum is 
observed, instead for a small excitation spot the conden- 
sation occurs on a ring in momentum space. This differ- 
ence has been explained within a mean field theory based 
on the Gross-Pitaevskii equation, including driving and 
dissipation [io| . More recently, another remarkable phe- 
nomenon related to the flow in a continuously pumped 
polariton condensate was observed experimentally 
vortices are spontaneously created in polariton conden- 
sates without setting the system into rotation. A theo- 
retical interpretation of this effect was given in the frame- 
work of the generalized Gross-Pitaevskii equation. In a 
significant fraction of random landscape realizations the 
polariton condensate contains a vortex. A related predic- 
tion was made by Keeling and Berloff: they found that 
a rotating vortex lattice can be spontaneously generated 
in a large regular trap [l2| . 



The above mentioned phenomena can be understood 
within a mean field theory, i.e. a theory where the quan- 
tum polariton field is replaced by a classical field. In 
this approximation, all information on the fluctuations is 
however lost. Since we deal with a two-dimensional sys- 
tem, the physics of fluctuations in polariton condensates 
is in analogy with equilibrium systems expected to be 
very rich [7|, Il3l . |14| and the question arises for example 
to what extent the physics related to the Berezinskii- 
Kosterlitz-Thouless survives the driving and dissipation 
of polariton condensates. 
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FIG. 1: Sketch of the separation of phase space in two regions: 
the lower polariton and the exciton reservoir. Scattering from 
the exciton reservoir into the lower polariton region replen- 
ishes the lower polaritons. 

Experimentally, the fluctuations of the polariton con- 
densates have been investigated under the form of the 
first and second order coherence functions. In the first 
order equal time spatial correlation function, long range 
correlations were observed above the stimulation thresh- 
old for condensation [3, flah Other correlation functions 
include the temporal first |5[ and second order coherence 

The semiclassical Boltzmann equation [18], Ha, (20] pro- 
vides a theoretical description of the first order spatial 
coherence, which is the Fourier transform of the momen- 
tum distribution. Including the details of the relaxation 
mechanisms, this formalism is expeced to give a reliable 
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estimate for the required polariton density to achieve con- 
densation. Above the condensation threshold, the ran- 
dom phase approximation contained in the Boltzmann 
equation breaks down and more sophisticated techniques 
should be used. Schemes that have been implemented 
in the literature involve the separation of the condensate 
mode from the excited states [2l[ and a generalization 
of the Boltzmann equation that includes the coherences 
within a Bogoliubov approximation [22j |. 

One of the remarkable consequences of the nonequi- 
librium nature of the polariton condensates is that the 
collective excitation spectrum is changed at low wave vec- 
tors: a diffusive instead of soundlike character is found 
for the low energy phase modes. This dispersion of ele- 
mentary excitations was found in a theoretical descrip- 
tion based on a Keldysh Green function technique worked 
out by the Littlewood group [23[ . The same spectrum is 
straightforwardly recovered by linearizing the generalized 
Gross-Pitaevskii equation around a steady state [24|, a 
calculation that is easily extended to spatially nonuni- 
form situations. 

It is well known from quantum optics and the theory of 
weakly interacting Bose gases that fluctuations can be in- 
cluded by introducing a stochastic element in the Gross- 
Pitaevskii equation 25, 26, 27|,l28|. For polariton systems 
in the parametric oscillation regime, such a method was 
used in Ref . [2^| . One of the great advantages of these so 
called classical field methods is that the nonuniformity of 
a system does not introduce any appreciable extra cost 
in their numerical implementation. A second advantage 
is that practical numerical calculations do not require a 
perturbative expansion around a condensed state and can 
even be applied to stud y p hysics related to the conden- 
sation phase transition [30 . l3l| . Finally, these methods 
can describe the evolution of the system in real time so 
that information on both the steady state and transients 
can be obtained. The latter can be of particular use to 
model experiments that are performed under pulsed ex- 
citation [2|. 

Due to the approximations involved in the classical 
field methods, they cannot describe the particles up to 
arbitrary large momenta, where quantum effects (most 
notably spontaneous scattering) are dominant [25j |. In 
this respect, polariton condensates are very well suited 
for a classical field description, because, as illustrated in 
Fig. [U the phase space can be naturally divided into two 
parts: i) a low energy polaritonic region with a small 
effective mass that shows quantum coherence above a 
certain threshold density and ii) a high energy excitonic 
'reservoir' region with a high effective mass that under 
typical experimental densities behaves as an incoherent 
classical gas. The role of the two subsystems is very dif- 
ferent: the polaritonic field is the quantity of experimen- 
tal interest because it is easily accessible in photolumi- 
nescence experiments and can be driven into the quan- 
tum degenerate regime. The role of the reservoir is to 
replenish the polariton region through relaxation. 

We will present in this paper a set of classical field 



equations for the (Wigner distribution function of the) 
polariton dynamics coupled to the exciton reservoir and 
apply them to calculate the equal time first and second 
order coherence functions across the condensation thresh- 
old. For simplicity, we have not included the polariza- 
tion degree of freedom. A Boltzmann description of po- 
lariton condensates including polarization can be found 
in Ref. [321] . At thermal equilibrium, the magnitude of 
fluctuations can be parametrized by a single quantity, 
the temperature. For a weakly interacting system, this 
temperature can be extracted by fitting the tail of the 
momentum distribution with a Maxwellian curve. Ex- 
periments have shown that the Maxwell-Boltzmann dis- 
tribution is also recovered for the tail of the out of equi- 
librium polariton distribution 0, 0, H[ and has even been 
experimentally observed in the case of weak coupling las- 
ing [33l ] . We will show that out of equilibrium, the uni- 
versal characterization of fluctuations by a temperature 
parameter breaks down. We will point out several cru- 
cial aspects of the condensate-reservoir interactions that 
affect the correlation functions without changing the tail 
of the momentum distribution. 

We start by presenting the model Hamiltonian for the 
nonresonantly excited polariton system in Sec. [TTJ It is 
shown in Sec. HHI how a master equation for the lower po- 
lariton field can be derived and how to solve it within the 
truncated Wigner approximation in Sec. HVl The reser- 
voir dynamics is discussed in Sec. [V] In Sec. IVI1 we dis- 
cuss the relation between our model and the Boltzmann 
equation. Numerical Monte Carlo results are presented 
in Sec. IVIII Conclusions are drawn in Sec. IVIII1 

II. HAMILTONIAN 

In order to treat these two regions of polaritonic phase 
space with a very different character, we replace the orig- 
inal Hamiltonian by a Hamiltonian for polaritons and 
excitons, that are anihilated by the operators "0( x ) an d 
<^(x) respectively. In terms of these annihilation opera- 
tors, our model Hamiltonian reads 

H = Jdx [H LP { X ) + H R { X ) + ff fl ,ip(x)] . (1) 

The lower polariton Hamiltonian density is the usual 

H LP (x) = ^t (x) _Z^( x ) + |v t (x)V t (x)V(x)^(x), 

(2) 

where m^p is the lower polariton effective mass and g 
quantifies the strength of the polariton-polariton inter- 
actions, that is well approximated by a zero-range po- 
tential. The exciton reservoir Hamiltonian is given by 

ff R (x) = 0t (x) zZ! 0(x) + !^t (x ^t( x )^ x) 0( x) , (3) 
2m x 2 

In the polariton/exciton basis, the exciton-exciton 
Coulomb scattering gives rise to various coupling terms. 
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The relevant ones are 

H R , LP (x) = <^p(x) + fff£p(x) + H£ LP (x). (4) 
Lower polaritons are created by the term 

^l a Lp( x )=5^(x)^(x)0(x)^( x)j ( 5 ) 

whereas they are destroyed by 

<!p(x) = 9 t (x)0t( x )0( x ) V ,(x). (6) 

Mean field shifts of the lower polaritons due to the 
excitons in the reservoir and vice versa, are described by 
the Hamiltonian 

^lp(x)=3^(x)^(x)^(x)^(x). (7) 

Note that we have in fact extended phase space by in- 
troducing two particles: the excitonic phase space is ex- 
tended down to k = and the polaritonic phase space to 
arbitrarily large momenta. Both extensions however add 
a very few states in the physically relevant regions. 

III. MASTER EQUATION 

In order to take advantage of the incoherent nature 
of the excitons in the reservoir, we will trace them out 
from the dynamics and obtain a quantum equation for 
the LP field alone. The LP field dynamics can be studied 
through the Liouville equation for the density matrix 
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-i[H,p]. 



(8) 



Going through the usual steps in the derivation of the 
Master equation in quantum optics, the Master equation 
for the full density matrix reads in the second Born ap- 
proximation in Hr lp 



P (t) = p(t ) - i / dt'[H RiL p(t'),p LP {t')] 

Jto 



f dt' f dt"[H R , LP (t% [H R , LP {t"),p(t")]l (9) 

JtQ JtQ 



where H R ^p(t) is in the interaction picture with respect 
to the Hamiltonian Hq = Hlp + Hp. 

Taking the trace of this equation over the reservoir 
degrees of freedom gives the desired master equation for 
the reduced density matrix p^p of the lower polariton 
subsystem. The second term on the RHS vanishes when 
the trace over the reservoir is taken, so we only have to 
analyze the third one. It consists of terms like 



;ain/loss Tj-gain/loss 
R,LP R.LP 



9}, 



(10) 



where 



^ R {HTlpHTlp9lp} = Tr R {H^lpH^lpp LP } = 0. 

(11) 



A nonzero term is e.g. given by 

Ri = J dt' dt"Tr R {H R %(t')H^l P (t") p(t")}. (12) 

In order to work out the trace over the reservoir, we 
introduce relative and center of mass coordinates 



X = 

T = 



2 ' 
t' + t" 



t = t' - f' 



(13) 
(14) 



We define the Wigner transform of the reservoir propa- 
gator as 



„iujt „— ikx 



F (X, k, T,u) = J dt dx e luJt e 
Tr fl {0t(X + x/2, T + t/2)<j>(X - x/2, T - t/2)}. (15) 

With the inverse transformation, we obtain for Ri de- 
fined in Eq. ((I2|) 

El = ^3 E / dXdx e j(Ak x - Aet) n/CX.ki.a.a) 

kl.2,3 

^ (X + x/2, T + t/2)i/>(X - x/2, T - t/2)p(T - t/2), 

(16) 

where fi is the area of our system, Ak = k2 + k3 — ki , 
Ae = e(k 2 ) + e(k 3 ) — e(ki) and n/(X, ki j2 ,3) is a typical 
Boltzmann collision rate (density in phase space) 

n / (X,k lj3i3 ) = /(X,k 1 ,T)[/(X,k 2 ,T)+l][/(X,k 3 ,T)+l] 

_ (17) 

We have used the quasi-particle approximation [34j 

F W (X, k, T, w ) = (2Ki)8(u - 6 k )/(X, k, T). (18) 

The time evolution of the LP field operators is approxi- 
mately given by 

^(X+x/2,T+t/2) ~l^e iQ ! x+I / 2 )e^V(Q,J'-i/2). 



Q 



(19) 

where the interaction shift in the frequency of ipQ was ne- 
glected. The exponential e lC * x can be combined with the 
exponential in Eq. (fl6|) . Because the typical reservoir 
momentum is much larger than the typical lower polari- 
ton momentum (see Fig. [1}, this factor is negligible. For 
the same reason, also the x in the second field operator 
in Eq. ([To]) can be neglected. If we then also assume that 
the density matrix is slowly varying on the microscopic 
time scale t, the integral over the relative time imposes 
energy conservation for the scattering process. We can 
then finally rewrite Eq. (fT2|) as 



ki,k2,k3,Q 

n/(X, k h2 , 3 )^ Q (T)i> x (T)p(T). (20) 
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The main simplifying assumption of the model consists 
now of assuming that the expression ([20)) is a function of 
the total reservoir density n R and the energy e q only. 
This comes down to the assumption of a steady state 
distribution of the reservoir excitons among the different 
k states. 

Working out the trace over the reservoir in Eq. (J9j) 
yields gain and loss terms for the lower polariton field 
from the collisons involving reservoir excitons. Collecting 
all these terms, we obtain 

j t p(t) = -i[H LP ,p]+K in (p) + K out (p), (21) 

where the density matrix evolves under the in-scattering 
as 

Kin(p) = \Y,J d*Rrn{n R ,e n ) [ e "^(x)W(q) 

-e^(q)^(x)p + h.c] , (22) 
and under out-scattering as 

^out(p) = \Y,j ^out(«fl,e q ) [e^q^M 

-6 i ^( X )^(q)p + / l .c] , (23) 

The rates -Rin/out are given by the usual semiclassical 
Boltzmann rates. Neglecting stimulated processes in the 
reservoir, R ln and i? ut depend on the reservoir density 
respectively as n 2 R and n R . We therefore write 

Rin(n R ,eq) = n^-R in (e q ) (24) 
RovLt(nR, e q ) = n R R out (6q) (25) 

Actually, another loss mechanism for the lower polari- 
ton field is present: leakage of the photon out of the im- 
perfect microcavity mirrors, that gives a finite line width 
7 to the lower polariton. This loss mechanism has a neg- 
ligible energy and momentum dependence and can be 
added to the model by simply adding the constant term 
7 to R out . 



IV. WIGNER 

An exact solution of the Master Equation (f2Tj) is not 
possible, but numerical progress can be made by the use 
of quasi-probability distributions from quantum optics. 
In the presence of dissipation, the Wigner distribution 
function is believed to give robust results (see Ref. [35j . 
pps. 115,124). This method has been applied to study 
BEC aspects of parametrically generated signal polariton 
in microcavities in Ref. [29|. 

The Wigner distribution function is a quasi probability 
distribution defined on the space spanned by the complex 
valued functions "0(x). In order to avoid ambiguity, we 
will use below an explicit 'hat' notation for the quantum 
field operator V»(x). 

In terms of the density matrix, the Wigner distribution 
function is defined as 

JVhKx)] = ^ yd 2 A(x)exp[^(x)A(x)* -^(x)*A(x)] 

d 2 a(x)(a(x)| j oexp[A(x)^ t (x)-A*(x)v3(x)|a(x))], 

(26) 

where |a(x)) is a coherent state of polaritons at posi- 
tion x with complex amplitude a(x). Expectation values 
calculated with the Wigner distribution function corre- 
spond to expectation values of symmetrized operator ex- 
pressions. For example for the one-bbody density matrix, 
we have: 

J dty(x) PwM*)] r(x)v>e>o 

= -TrM^x^x') + ^( x ')V3t(x)}. (27) 

Using the operator correspondences [351 ] . the equation 
of motion for the Wigner quasi-probability distribution 
Pw is computed: 



aiy[y.(x),y.»(x)] 

dt 



d 



-F, 



d^y det d^*{*y det ' 9v(x)V'*(x) 



F, 



2 



-[7 + ft in +7e out ] 



+i 



if 



2AF dip(-x)ip*(x) 



d 



V*(X) - TT-T^X) 



JV[^(x),^*(x)]. (28) 



F^et is the deterministic mean field force acting on the polaritons 



Fdet = ~i 



h 2 V 2 : i{K in -K out -j) , g 



2m 



^(x). 



(29) 
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In Eq. (|28|) , a momentum cutoff for the field ip is im- 
plicitely introduced by formulating the problem on a spa- 
tial grid with cell area AT^. The expression 7^i n ,outV' 
should be understood as 

7e in V>(x) - r4(x)^ e -^ x 'i? in (e q )^xs (30) 
q 

HtoVM = n R (x)^e- icix ' J R out (e q )Vx'- (31) 
q 

From the mathematical point of view, the last term in the 
equation of motion (|28j) has proved to be an insurmount- 
able problem. If this term is neglected (the so-called 
truncated Wigner approximation), the quasi-probability 
distribution Pyy obeys a standard Fokker-Plank equa- 
tion, that correponds to the Langevin equation 

tty(x) = F det [-0(x),-0*(x)] +dW(x), (32) 

where dW is a complex Gaussian stochastic variable with 
the correlation functions: 



(dW(x)dW{x)) = 0, 
dt 



(dW(x)dW*(x')) = 



AAV 



((x|ftf n + 7^ ut |x'} + 2 7( 5 x , x /), 



(33) 



where U^ out = [7£ in ,out + (^in,out) T ]/2 are the sym- 
metrized kernels. 

Let us now estimate the order of magnitude of the 
third order derivative in Eq. (j2"5j) with respect to the 
other terms, in particular the second order derivative 
terms. The function P\y is peaked around the value of 
the field ip(x) whose squared modulus equals |?A(x)| 2 = 
iV(x) + 1/2. The variation of Pw occurs on a scale of 
its argument of order one. Derivatives are therefore ex- 
pected to be of the order of the function Pw itself and 
the prefactors determine the relative importance of the 
derivative terms in Eq.(4). This leads us to the con- 
clusion that the third order derivative is negligible with 
respect to the second order one if 



7 > 
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AV 



(34) 



The dissipative character of the system thus increases the 
region of validity of the truncated Wigner approximation. 
The dissipation gives away information about the system 
and destroys nontrivial quantum states (e.g. number or 
Schrodinger cat states). In terms of the Wigner function, 
oscillations of Pw accompanied by regions where it be- 
comes negative (that cannot be represented by a regular 
probability distribution) are washed out by the dissipa- 
tion (Hi. 



V. THE EXCITON RESERVOIR 

In our description of the microcavity dynamics, the 
exciton-like particles arc treated as a classical reservoir. 



This approximation allowed to trace out the excitonic 
degrees of freedom and to isolate the quantum dynamics 
of the polaritons from the classical exciton dynamics. In 
principle, the reservoir density appears as a deterministic 
classical quantity in the resulting equations of motion 
for the lower polariton dynamics. Physically, however, 
this is not expected to be a very good approximation, 
because the condensate serves as a relaxation mechanism. 
Stimulated scattering makes the rate of this relaxation to 
depend on the reservoir population. Similar ideas have 
been implemented in Ref. [2l| . where the dynamics of 
a single condensate mode was coupled to a Boltzmann 
equation for the excited states, and in Ref. [5(, where the 
reservoir was modelled by a saturable gain medium, a 
model widely used in laser physics (4l| . 

We propose to go beyond the approximation that the 
reservoir is unaltered by the system by coupling its dy- 
namics to the equation of motion for the classical polari- 
ton field 



dn R 
dt 



-ln[nR - n° R (I p ,ip)], 



(35) 



where n° R (I pi ip) is the average steady state value of the 
reservoir density in the presence of a pump with intensity 
I p and a lower polariton field ip. The relaxation time 
7^ is a measure of the time it takes for the reservoir 
density to adjust to a new environment (I p ,ip). Spatial 
diffusion of the reservoir excitons is expected to be a small 
effect [24[ and was therefore neglected. For the steady 
state value of the reservoir density, we assume that it 
is simply proportional to the balance of incoming and 
outgoing particles 



n° R (P,i>)=P(I p --&H)\ ies ) 



(36) 



where ^(VV)lrea = 2Re[V>*(ft in - TZoutM is the net 
scattering rate from the reservoir into the lower polariton 
branch. It is instructive to substitute Eq. (f36|) into Eq. 



^^P-lBUR-hRj^H)^ (37) 

where P = f3I p is the effective pump term for the ac- 
tive reservoir polaritons. The parameter /3 quantifies the 
backaction of the condensate on the reservoir. This back- 
action is needed to obtain a steady state for the dynami- 
cal equations above the threshold, where for ur = P/~fR 
the in-scattering rate exceeds the out-scattering rate. In 
mean field theory, the reservoir density ur is clamped to 
its threshold value «_R,m/ that satisfies for homogeneous 
systems n 2 Rmf R in (0) - n Rim fR out (0) = 7. If we rewrite 
the motion equations for tir in terms of the renormalized 
n R = n R /n R: , nf , we have 



dfiR 
~dT 



= P - iRflR 



(38) 



where a = f3jR/nR^ m f. Also in the presence of fluctua- 
tions, the dimensionless reservoir density fiR is close to 
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one above threshold, in order for the gain to compensate 
for the losses. The factor a plays an important physi- 
cal role because the backaction of the condensate on the 
reservoir tends to damp the condensate fluctuations. If 
the condensate density is at some time larger than aver- 
age, the reservoir will be depleted, i?; n — -Rout decreases 
and the deterministic part in the equations of motion for 
the condensate will decrease the amplitude of the fluc- 
tuation. In principle, the parameter a could be calcu- 
lated from the Boltzmann equation. We prefer however 
to study the physics in terms of this parameter, because 
it gives a good insight in the nonequilibrium aspects of 
the coherence. 

In the truncated Wigner approximation, the density 
of polaritons is related to ip as n — \tp\ 2 — 1/(2AV), 
or in words, the classical field ip contains half a particle 
per mode of zero point fluctuations. These fluctuations 
should be taken into account when evaluating the last 
term in Eq. (J38J. For the out-scattering, the zero-point 
fluctuations do not contribute and should be subtracted, 
whereas for the in-scattering, the zero-point fluctuations 
give rise to only half of the spontaneous in-scattering. 
The remaining part should be added. The equation of 
motion for the reduced reservoir density then finally reads 

=P~ IR^R - a 4; (V> f 1p) \res,W 

at at 



2AV 



^[i?out(e k )+i?out(ek)]. (39) 



VI. 



RELATION WITH THE BOLTZMANN 
EQUATION 



In the dissipative case, the derivation of the truncated 
Wigner equation did not rely on the formation of a con- 
densate. We can therefore describe within the same for- 
malism the condensed and non-condensed polariton gas. 
In the case that coherence is negligible, it is instructive 
to simplify the stochastic equations of motion ([32]) . We 
will find that in the incoherent regime, the polariton con- 
densate can be described with a Boltzmann-like equa- 
tion d3. 

For simplicity, we consider the case of a uniform reser- 
voir density. By writing the stochastic motion equations 
for the field ipQs) in momentum space, treating the inter- 
actions in the second Born approximation, and assuming 
that there are no phase relations between the different 
momentum components {tp*(k, t)ip(k! , t)) — [N(k,t) + 
l/2]<$k.k'; one obtains the following Boltzmann like equa- 
tion for the time evolution of the densities in momentum 
space 



dN(k, t) 
dt 



= flin(ek) [N(k, t) + 1] - [i? out (e k ) + j]N{k, t) 
+ I B [N(k)]+I Q [N(k)}. (40) 



losses through the cavity mirrors. Collisions are de- 
scribed by the last two terms. Ib is the usual Boltzmann 
collision integral: 

I B [N(k)} = -2ng 2 £ S^ R (e 1 + e 2 - e 3 - e 4 ) 

ki,k 2 

x [NxN 2 (l + 2V 3 )(1 + N 4 ) - N 3 N 4 (l + 2Vi)(l + 2V 2 )] , 

(41) 

wheveN 1 (t) = N(k,t),N 2 (t) = N(k 1 +k 2 -k,t),N 3 (t) = 
-ZV(ki, t), N^t) = N(k 2 ,t) and analogous for the en- 
ergies ti. The 8- function for energy conservation is 
broadened due to the finite lifetime of the polaritons: 
5 v (lo) = s\a{ujv) / (ttlo) . 

The extra collisional term is due to to the fact that 
our stochastic classical field model does not coincide with 
the true quantum dynamics (the third order derivatives 
in Eq. (j2"5)) are neglected): 

2 

Ic[N(k)} = *7-«(ci + £ 2 ~ £3 - £4) 

ki,k 2 

x [JVi+JVb-JVg-jvy. (42) 

This term is spurious, because the Boltzmann equation 
should be recovered in the incoherent limit. Our classical 
field model can therefore only be a good approximation 
of the full quantum dynamics if the term Ic is negligible 
with respect to the other terms in Eq. (|4TH) . It scales 
as Ic oc (g/ AV)(gnAV). If the occupation numbers per 
grid cell nAV are much larger than unity, the Boltzmann 
collision term Ig is obviously dominant with respect to 
Ic ■ This is the typical condition for the use of the Wigner 
distribution function in the description of a stable Bose 
gas [2^. For bosons with a finite life time, this condition 
can be relaxed, because even when nAV is not much 
larger than unity, the spurious term can be still much 
smaller than the reservoir term i?; n . The in-scattering 
rate Ri n should compensate the losses 7. If occupation 
numbers are not large, the truncated Wigner is therefore 
still expected to yield physical results if g/AV <C 7. Note 
that the latter requirement coincides with the condition 
(f34|) derived from the full equation of motion (1281) . 

If we neglect the collisional terms in Eq. (|40|) . the 
steady state solution is 



N(k) = 



1 



[7 + #out(ek)]/#m(ek)-l' 



(43) 



The simplest model that yields a temperature Tr (that is 
in experiments typically higher than the lattice tempera- 
ture) for the tail of the polariton momentum distribution 
is obtained by setting the out-scattering rate to zero 



Rout(E) — 0, 
R in (E) oc 7 exp(- J E;/T H ). 



(44) 
(45) 



The first two terms describe the evolution of the mode 
occupation due to the interaction with the reservoir and 



Studies of the Boltzmann equation [381 ] have however 
shown that the rates Ri n and i? ou t both tend to incre- 
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asse as a function of the energy, approximately as 

i? in oc exp{E/k B T R ), (46) 
R out tx exp(2£/fe B T R ). (47) 

We will see below that a nonzero out-scattering enhances 
fluctuations. 



VII. NUMERICAL RESULTS 

The stochastic motion equations can be simulated by 
Monte Carlo techniques. As we have already mentioned 
in the introduction, the nonequilibrium condition of the 
polariton condensates makes that the effect of the reser- 
voir on the correlation functions cannot be characterized 
by the temperature alone. We will discuss below two 
other physical quantities that determine the degree of 
coherence in the polariton condensate: the feedback pa- 
rameter a and the out-scattering rate i? ut- The other 
parameters, we keep fixed for all simulations: m/h = 
l^m -2 meV _x , g/h = 0.03^m 2 and k B T R = 2meV. The 
simulations were done on a 32 x 32 point grid with physi- 
cal dimension of 66 x 66/im 2 and periodic boundary con- 
ditions. 

Figs. [H and [3J illustrate single Monte-Carlo realiza- 
tions of the classical field ^(x). Even though these im- 
ages have strictly speaking no direct physical meaning, 
they already illustrate qualitatively the coherence prop- 
erties of the polariton condensate. Fig. [5] shows two 
examples for a finite excitation spot for pump intensities 
below (panels a and b) and above the threshold (panels 
c and d). At low density, both the density and phase 
fluctuations are large, whereas the phase fluctuations are 
clearly suppressed in the high density regime. Panel (d) 
shows that phase coherence exists all over the extent of 
the excitation region. The concentric phase profile orig- 
inates from the repulsive polariton-polariton interaction 
that causes an outward flow of polaritons [lOj . 

Figure [3] shows snap shots of the polariton density and 
phase for a uniform pump below and slightly above the 
threshold. The phase profile of panel (d) shows that the 
phase ordering is only partial. Vortex-anti vortex pairs 
appear to exist at densities well above the stimulation 
threshold. This is an indication that the physics of the 
Berezinskii-Kosterlitz-Thouless type could occur in po- 
lariton condensates. 

Three momentum distributions for increasing pump in- 
tensity are shown in Fig. 0J As expected, our model 
shows the build up of a large occupation of the low mo- 
mentum states for increasing pump intensity. The mo- 
mentum distributions appear to be rather well fitted by 
a Bose-Einstein function (full line). It is important to 
mention here the important role of the reservoir relax- 
ation rates. We have chosen them in such a way that a 
thermal distribution is obtained even in the absence of 
collisions between lower polaritons. In simulations with 
energy independent relaxation rates (not shown) and a 
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FIG. 2: Snapshots of a single Monte Carlo realization of the 
density (upper panels) and phase (lower panels) for a finite 
size excitation spot with intensity below (left hand panels) 
and above threshold (right hand panels). 
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FIG. 3: Snapshots of a single Monte Carlo realization of the 
density (a,c) and phase (b,d) for excitation parameters below 
(a,b) and above threshold (c,d). 



large, yet realistic [42j polariton-polariton interaction 
strenth, we have obtained a constant instead of exponen- 
tial decay at large momenta. 

Note that the temperature extracted from the fits of 
the tails to a Maxwellian is lower than the reservoir tem- 
perature Tr (2 meV for the present simulations), that 
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enters the rates i?i n ,out according to Eq. (|47j) : the non- 
linearity modifies the temperature that is expected in the 
linear regime. We remind the reader that Tr does not 
coincide with the lattice temperature and that Tgt < Tr 
does not imply that the polariton temperature is lower 
than the lattice temperature. 
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FIG. 5: Spatial coherence (open circles, the grey band indi- 
cates the error on the Monte Carlo data) for two values of the 
feedback from the condensate on the reservoir: from left to 
right a — 0.01/im 2 (a,c) and a = 0.1/im 2 (b,d) . The full line 
shows the decay of the correlations in a non-interacting Bose 
gas, at the temperature reported in the panels. 
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FIG. 4: Momentum distribution of the polariton field for var- 
ious pump intensities. The open circles are the result of the 
Monte Carlo simulations and the full line is a fit to a BE dis- 
tribution. Statistical errors of the Monte Carlo simulations 
are within the symbol size. 

The subtle features of long range coherence are much 
clearer in the Fourier transform of the momentum distri- 
bution, i.e. the first order spatial coherence function. In 
Fig.[5]two values of a are compared. Below the condensa- 
tion threshold, the gain saturation parameter a (see Eq. 
(155)0 has no influence and the fit of the coherence by 
the of the noninteracting Bose gas yields the reser- 
voir temperature of 2 meV. For the simulation above the 
threshold, a higher value of a improves the long range co- 
herence. Both spatial coherence functions are relatively 
well fitted by the one of a noninteracting Bose gas. Both 
temperatures are below the reservoir temperature. The 
lowest effective temperature is obtained for the largest 
feedback parameter a. 

In the simulations of Fig. the out-scattering was 
set to zero. In the simulations presented in Fig. [6l we 
have included this effect. In order to avoid exceedingly 
large rates in the model, we have put a cutoff in the 
magnitude of R out as R out {E) = mm[e 2E ^ Ta , 3.3meV]. 
The in-scattering rate was chosen Ri n (E) = [R ou t(E) + 

Fig. [5] shows that the out-scattering has a big effect 
on the coherence function. This should not come as a 



surprise, because the out-scattering increases the fluctu- 
ations (physically shot noise due to the discrete nature of 
the polariton field). Keeping a = O.Ol^m -1 as in Fig.0 
but including some out-scattering, the coherence in panel 
(a) is dramatically decreased. As compared to the simu- 
lations of Fig. the effect of a is much more pronounced. 
For the smallest value of a, the temperature of 5.5 meV 
for which the spatial coherence is reasonably well fitted, 
is much larger than the one that is extracted from the tail 
of the momentum distribution (less than 1 meV): the po- 
lariton condensate behaves in this regime very different 
from the ideal Bose gas. 



k (a) 






T=5.5 meV 



10 20 30 10 20 30 

x(^im) x((im) 

FIG. 6: The spatial coherence function as in Fig. [5] but in- 
cluding out-scattering. Simulatios with a — 0.01/zm 2 (a) and 
a = 0.1/rni 2 (b) are shown. 

Another quantity of great physical interest is the 
second order coherence function (x, t; x', t') = 
(-0 f (x, t)^(x',t')ip(x',t')ij(yi, t)), that quantifies the 
density fluctuations. Experimentally, the equal position 
second order coherence was investigated, in Ref. [l7j for 
equal times t — t' and in Ref. [5] as a function of the 
delay t — t'. Within the Wigner formalism, different time 
correlation functions are not straighforwardly calculable, 
so we present here only results for the equal time second 
order coherence. 

Results of the equal position second order coherence 
g( 2 )(0) = g( 2 )(x, i;x, t) are shown in Fig. [71 for several 
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parameter values. As expected, g^(0) approaches the 
value 2 of the incoherent Bose gas in the low density 
regime. For increasing polariton densities, the second 
order coherence decreases, but its actual value depends 
again strongly on the chosen parameter values. A larger 
value of the feedback parameter a, suppresses the den- 
sity fluctuations. This is in agreement with the model 
described in Ref. where the density fluctuations are 
proportional to the saturation density (large saturation 
density means small feedback from the condensate on 
the reservoir). Fig. [7] also shows that the out-scattering 
increases the density fluctuations. This dependence is 
expected, because adding the knock out processes leaves 
the deterministic term in the evolution equation for the 
classical field unaltered, but increases the fluctuations. 

Note that the density fluctuations are within our model 
not always monotonous, but for some parameter values 
show a minimum value slightly above threshold. Non- 
monotonous behavior of g^ 2 '(0) was also observed in ex- 
periments on polariton condensation in CdTe microcav- 
ities [13] ■ Also in the theoretical work of Ref. [H, H(| 
based on a Boltzmann equation for the excited states cou- 
pled to a master equation for the condensate mode, an 
increase of density fluctuations above the threshold was 
found. It is however important to mention that in Fig. [7] 
the interaction energy is very large when t/ 2 -* (0) increases 
again (1 meV blue shift due to condensate-condensate 
interactions alone). When the value of the blue shift is 
reduced to 0.2 meV, <?^(0) is found to be very close to 
one. 

We want to point out that we have not found any 
regime with good long range spatial coherence and large 
density fluctuations. Indeed, Fig. [5] (a) shows that at 
the density n w 20/im~ 2 the spatial coherence is, al- 
though longer than the thermal de Broglie wave lenth 
corresponding to Tr, limited to about 10 /im. Physically 
it is actually not expected that good spatial coherence 
and strong density fluctuations can go together, because 
phase fluctuations are coupled to the density fluctuations 
through the interaction and kinetic energy. So far, in the 
experiments on CdTe microcavities where the increase of 
(0) as a function of pump power was observed, no de- 
crease in spatial coherence was seen. It is possible that 
the distance at which the spatial coherence was probed 
is too short for the decrease in spatial coherence to be 
detectable, but we cannot exclude other explanations in 
terms of extrinsic experimental effects. The measured 
density fluctuations could for example contain a compo- 
nent due to intensity fluctuations in the excitation laser. 

VIII. CONCLUSIONS 

We have derived classical field equations for a nonres- 
onantly excited polariton condensate in a semiconductor 



microcavity within the truncated Wigner approximation. 
Thanks to the polariton losses our model remains physi- 
cal in the low density regime and allows to describe the 
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FIG. 7: Second order coherence function as a function of the 
total density obtained for different excitation powers and the 
simulation parameters indicated in the legend. 



polariton condensate at all densities. Our equations were 
shown to reduce to the Boltzmann equation in the low 
density regime below threshold. Above threshold, the 
equations were analyzed numerically with Monte Carlo 
simulations. The first and second order spatial coherence 
were shown to depend dramatically on the feedback from 
the condensate on the reservoir (the gain saturation) and 
on the collisions with reservoir excitons that knock po- 
laritons out of the condensate. Within our model, the 
density fluctuations can show nonmonotonous behavior 
as a function of the polariton density. We predict that 
an increase in density fluctuations is accompanied by a 
decrease in the spatial coherence. 

Finally, the vortex defects in individual Monte Carlo 
realizations of the polariton field show that the spatial 
coherence is limited by the spontaneous appearance of 
vortex defects in the phase. A further study of the role 
of vortices is necessary to understand their effect on the 
spatial coherence. 
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